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Preface 



Hamiltonian Boundary Value Methods (in short, HBVMs) is a new class of nu- 
merical methods for the efficient numerical solution of canonical Hamiltonian 
systems. In particular, their main feature is that of exactly preserving, for the 
numerical solution, the value of the Hamiltonian function, when the latter is a 
polynomial of arbitrarily high degree. 

Clearly, this fact implies a practical conservation of any analytical Hamilto- 
nian function. 

In this notes, we collect the introductory material on HBVMs contained in the 
HBVMs Homepage, available at the url: 

http : / / web .matji . unif i . it /users /brugnano/HBVM/ index . html 
The notes are organized as follows: 

• Chapter 1 : Basic Facts about HBVMs 

• Chapter 2: Numerical Tests 

• Chapter 3: Infinity HBVMs 

• Chapter 4: Isospectral Property of HBVMs and their connections with Runge- 
Kutta collocation methods 

• Chapter 5: Blended HBVMs 

• Chapter 6: Notes and References 

• Bibliography 
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Chapter 1 

Basic Facts about HBVMs 



We consider Hamiltonian problems in the form 

y{t) = JVH{y{t)), y{to) = yo e M^'-, (1.1) 

where J is a skew-symmetric constant matrix, and the Hamiltonian H{y) is as- 
sumed to be sufficiently differentiable. Usually, 

so that (11.11) assumes the form 

q = VpH{q,p), p = -VgH{q,p). 

The induced dynamical system is characterized by the presence of invariants of 
motion, among which the Hamiltonian itself: 

H{yit)) = VH{y{t)Yy{t) = V/J(y(t))^JV//(y(t)) = 0, 

due to the fact that J is skew- symmetric. Such property is usually lost, when 
numerically solving problem (11.11) . This drawback can be overcome by using 
Hamiltonian BVMs (hereafter, HBVMs). 

The key formula which HBVMs rely on, is the line integral and the related 
property of conservative vector fields: 

H{y^) - H{yo) = h [ a{t^ + Thf\/H{a{to + r/i))dr, (1.2) 

for any yi E M^"^, where a is any smooth function such that 

o-(^o)=2/o, a{to + h)=yi. (1.3) 



5 



6 
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Here we consider the case where a(t) is a polynomial of degree s, yielding an 
approximation to the true solution y(t) in the time interval [to, to + h]. The nu- 
merical approximation for the subsequent time-step, yi, is then defined by (11.31) . 
After introducing a set of s distinct abscissae 

< ci,...,c, < 1, (1.4) 

we set 

Yi = o-{to + Cih), i = l,...,s, (1.5) 

so that cr(t) may be thought of as an interpolation polynomial, interpolating the 
fundamental stages Yi, i = 1, . . . , s. We observe that, due to (|1.3I) . a{t) also 
interpolates the initial condition yo. 

Remark 1. Sometimes, the interpolation at to is explicitly required. In such a 
case, the extra abscissa Cq = is formally added to This is the case, for 

example, of a Lobatto distribution of the abscissae 

Let us consider the following expansions of a{t) and a{t) for t E [to, to + 

cr(to + rh) = V^Tj^il^), (^{to + ^h) = l/o + ^ Y^Tj / Pji^) dx, (1.6) 

where {Pj{t)} is a suitable basis of the vector space of polynomials of degree at 
most s — 1 and the (vector) coefficients {7j} are to be determined. Because of the 
arguments in [61171 [H, we shall consider an orthonormal basis of polynomials on 
the interval [0, 1], i.e.: 

/ P;{t)P,{t)dt = 8,„ i,j = l,...,s, (1.7) 
Jo 

where 5ij is the Kronecker symbol, and Pi{t) has degree i — 1. Such a basis can 
be readily obtained as 



/^(t) = v/2r^P,_i(t), i = l,...,s, (1.8) 

with Pi-i{t) the shifted Legendre polynomial, of degree z — 1, on the interval 
[0,1]. 

Remark 2. From the properties of shifted Legendre polynomials (see, e.g., [T] or 
the Appendix in /jdj/j, one readily obtains that the polynomials {Pj{t)} satisfy the 
three-terms recurrence: 

P,{t) = 1, P2(t) = \/3(2t-l), 
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We shall also assume that H{ii) is a polynomial, which implies that the inte- 
grand in (11.21) is also a polynomial so that the line integral can be exactly computed 
by means of a suitable quadrature formula. In general, however, due to the high 
degree of the integrand function, such quadrature formula cannot be solely based 
upon the available abscissae {cj}: one needs to introduce an additional set of ab- 
scissae {ci, . . . , Cj.}, distinct from the nodes {cj}, in order to make the quadrature 
formula exact: 

f &{to + ThfVH{a{to + Th))dT= (1.9) 
Jo 

s r 

J2 Pi^to + CihfVH{(x{to + Cih)) + Pi&{to + CihfVH{a{to + q/i)), 

i=l 1=1 

where i = 1, . . . , s, and i = 1, . . . , r, denote the weights of the quadrature 
formula corresponding to the abscissae {c,} and {cj}, respectively, i.e., 

(1.10) 

Remark 3. In the case considered in the previous Remark [7] i. e. when cq = 
is formally considered together with the abscissae di.4D . the first product in each 
formula in U.IO^ ranges from j = to s. Moreover, also the range of {(3i} 
becomes i = 0, 1, . . . , s. However, for sake of simplicity, we shall not consider 
this case further 

According to [|28l . the right-hand side of (11.91) is called discrete line integral, 
while the vectors 

Yi = a{to + c,h), i = l,...,r, (1.11) 

are called silent stages: they just serve to increase, as much as one likes, the 
degree of precision of the quadrature formula, but they are not to be regarded as 
unknowns since, from (|1.6I) . they can be expressed in terms of linear combinations 
of the fundamental stages (11.51) . 

Definition 1. The method defined by substituting the quantities in (11.61) into the 
right-hand side of (11.91) . and by choosing the unknown coefficients {7^} in or- 
der that the resulting expression vanishes, is called Hamiltonian Boundary Value 
Method with k steps and degree s, in short HBVM(A;,s), where k = s + r 
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In such a way, one easily obtains, from (ll.2l) - (ll.3l) . 

H{a{to + h)) = H{yo), 

that is, the value of the Hamiltonian is exactly preserved at the subsequent approx- 
imation, provided by cr(to + h). 

In the sequel, we shall see that HBVMs may be expressed through different, 
though equivalent, formulations: some of them can be directly implemented in a 
computer program, the others being of more theoretical interest. 

Because of the equality (11.91) . we can apply the procedure directly to the origi- 
nal line integral appearing in the left-hand side. With this premise, by considering 
the first expansion in (11.61) . the conservation property reads 

VtJ / P,{r)VH{a{to + Th))dT = 0, (1.12) 
i=i 

which, as is easily checked, is certainly satisfied if we impose the following set of 
orthogonality conditions 

^•=/ PAr)JVH{a{to + Th))dT, j = l,...,s. (1.13) 
Jo 

Then, from the second relation of (11.61) we obtain, by introducing the operator 

L{f;h)a{to + ch)= (1.14) 

PC fl 

a(to) + /iV/ Pjix)dx / P,(r)/((T(to + r/i))dr, ce[0,l], 
j^i Jo Jo 

that a is the eigenfunction of L{JVH; h) relative to the eigenvalue A = 1: 

a = L{JVH]h)a. (1.15) 

Definition 2. Equation ( 17.751) is the Master Functional Equation defining a 

Remark 4. From the previous arguments, one readily obtains that the Master 
Functional Equation ( 17.751) characterizes HBVM{k, s) methods, for all k > 1. 
Indeed, such methods are uniquely defined by the polynomial a, of degree s, the 
number of steps k being only required to obtain an exact quadrature formula (see 

To practically compute a, we set (see (11.51) and (11.61) ) 

s 

Yi = a{to + Cih) = yo + h^aij'jj, i = l,...,s, (1.16) 



9 



where 



^ij= / Pj{.x)<lx, i,j = l,...,s. (1.17) 

Inserting (11.131) into (11.161) yields the final formulae which define the HBVMs 
class based upon the orthonormal basis {Pj}'- 



s „i 

Yi = yQ + hy^ai, / P,(r)JV/f(a(to + r/i))dr, 

.7=1 



i = l,...,s. (1.18) 



For sake of completeness, we report the nonlinear system associated with the 
HBVM(A;, s) method, in terms of the fundamental stages {Yi} and the silent stages 
[Yi] (see (|l.lll) ). by using the notation 

f{y) = JVH{y). (1.19) 

In this context, it represents the discrete counterpart of (11.181) . and may be directly 
retrieved by evaluating, for example, the integrals in (11.181) by means of the (exact) 
quadrature formula introduced in (11.91) : 

Y, = (1.20) 

s / s r \ 

3=1 \l=l 1=1 J 

From the above discussion it is clear that, in the non-polynomial case, supposing 
to choose the abscissae {q} so that the sums in (11.201) converge to an integral as 
r = /c — s — )■ oo, the resulting formula is (|1.18l) . This implies that HBVMs may 
be as well applied in the non-polynomial case since, in finite precision arithmetic, 
HBVMs are indistinguishable from their limit formulae (11.181) . when a sufficient 
number of silent stages is introduced. The aspect of having a practical exact 
integral, for k large enough, was already stressed in [l3l 161171 12411281 . 

We emphasize that, in the non-polynomial case, (11.181) becomes an operative 
method, only after that a suitable strategy to approximate the integral is taken 
into account. In the present case, if one discretizes the Master Functional Equa- 
tion (I1.14I) - (I1.15I) . HBVM(A;, s) are then obtained, essentially by extending the 
discrete problem (11.201) also to the silent stages (11.1 II) . In order to simplify the 
exposition, we shall use (11.191) and introduce the following notation: 

{r,} = {q} U {q}, {uJi\ = {A} U 

(1.21) 

yi = cr{to + Tih), fi = f{a(to + Tih)), i = l,...,k. 
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The discrete problem defining the HBVM(A;, s) then becomes, 

yi = yo + h^ Pj{x)dx^uJePj{Te)fi, i = l,...,k. (1.22) 

i=l £=1 

Remark 5. We also observe that, from l\1.13\) and the first relation in 47 .(51) . one 

obtains the equations 

s „i 

(T(to + nh) = V Pjin) / Pj{r)JVHiaito + r/i))dr, i = l,...,k, 



(1.23) 

which may be viewed as extended collocation conditions according to / [28l Sec- 
tion!], where the integrals are (exactly) replaced by discrete sums. 

By introducing the vectors 

?/ = (yf,...,y,T, e=(l,...,l)^GM^ 

and the matrices 

= diag(a;i, . . . , u;^), X„ e M'=^^ (1.24) 
whose (i, j)th entry are given by 

(X,),, = r P,{x)dx, (Vsh, = P,ir,), (1.25) 

we can cast the set of equations (|1.22l) in vector form as 

y = e®yo + h{IsVjn) ® Is™ f{y), (1.26) 

with an obvious meaning of f{y). Consequently, the method can be seen as a 
Runge-Kutta method with the following Butcher tableau: 



(1.27) 

















Ui ... Uk 



Remark 6. We observe that, because of the use of an orthonormal basis, the role 
of the abscissae {cj} and of the silent abscissae {cj} is interchangeable, within the 
set {tj}. This is due to the fact that all the matrices Is, Vs, and depend on all 
the abscissae {ti}, and not on a subset of them and, moreover, they are invariant 
with respect to the choice of the fundamental abscissae {q}. 
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The following result then holds true. 

Theorem 1. Provided that the quadrature defined by the weights {ui} has order 
at least 2s (i.e., it is exact for polynomials of degree at least 2s — 1), HBVM(k,s) 
has order p = 2s = 2 deg(cr), whatever the choice of the abscissae Ci, . . . , c^. 

Proof From the classical result of Butcher (see, e.g., [[221 Theorem 7.4]), 
the thesis follows if the usual simplifying assumptions C{s), B(p), p > 2s, and 
D{s — 1) are satisfied for the Runge-Kutta method defined by the tableau ([1.271) . 
By looking at the method ([1.26[) - ([1.27[) . one has that the first two (i.e., C{s) and 
B(p), p > 2s) are obviously fulfilled: the former by the definition of the method, 
the second by hypothesis. The proof is then completed, if we prove -D(s — 1). Such 
condition can be cast in matrix form, by introducing the vector e = G 
R*^\ and the matrices 

Q = diag(l, ...,s-l), D = diag(ri, . . . , r^), V = (Tf-^) G R^^^-\ 
(see also ([1.25[) ) as 

QV^Vt {isVjn) = (ee^ - v'^d) 

i.e., 

VslJnVQ = (ee^ - DV) . (1.28) 
Since the quadrature is exact for polynomials of degree 2s — 1, one has 

{iJnVQ)^. = |^5^c^,pP.(x)da;(jrr^)j = (^^' ^V,(x)dx(jt^-^)dt 

^ii- j Pi{x)x^dx^ , i = l,...,s, j = l,...,s-l, 

where the last equality is obtained by integrating by parts, with Sn the Kronecker 
symbol. Consequently, 



= (l-r/), i = l,...,fc, j = l,...,s-l, 
that is, ([1.28[) . where the last equality follows from the fact that 

J2 Piir) / Piix)x^dx = t\ j = 1, . . . , s - 1. □ 

£=1 
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Concerning the stability of the methods, the following result holds true. 

Theorem 2. For all k such that the quadrature formula has order at least 2s = 
2 deg{a), HBVM(k,s) is perfectly A- stable^ whatever the choice of the abscissae 

Ci, . . . , Cg. 

Proof As it has been previously observed, a HBVM(A;, s) is fully charac- 
terized by the corresponding polynomial cr which, for k sufficiently large (i.e., 
assuming that (11.91 ) holds true), satisfies the Master Functional Equation (11.141) - 
(11.151) . which is independent of the choice of the nodes ci, . . . , Cg (since we con- 
sider an orthonormal basis). When, in place of f{y) = JVH{y) we put the 
test equation f{y) = \y, we have that the collocation polynomial of the Gauss- 
Legendre method of order 2s, say cr^, satisfies the Master Functional Equation, 
since the integrands appearing in it are polynomials of degree at most 2s — 1, so 
that a = as- The proof completes by considering that Gauss-Legendre methods 
are perfectly A- stable. □ 

Example 1. As an example, for the methods studied in [^], based on a Lo- 
batto distribution of the nodes {cq = 0, ci, . . . , c^} U {ci, . . . , Ck-s], one has that 
deg((T) = s, so that the order of HBVM(k,s) turns out to be 2s, with a quadra- 
ture satisfying B{2k). Finally, we observe that, with such choice of the abscissae 
HBVM{s, s) reduces to the Lobatto IIIA method of order 2s. 

Example 2. For the same reason, when one considers a Gauss distribution for 
the abscissae {ci, . . . , Cg} U {ci, . . . , Ck-s}, cis done in one also obtains a 
method of order 2s with a quadrature satisfying B{2k). Similarly as in the pre- 
vious example, HBVM{s, s) now reduces to the Gauss-Legendre method of order 
2s. 

Remark 7. A number of remarks are in order, to emphasize relevant features of 
HBVM{k,s): 

• From Remark^ HBVM(k,s) are symmetric methods according to the time 
reversal symmetry condition defined in l \17i p. 218] (see also /fJ9l/). pro- 
vided that the abscissae {rj} (see A1.21\) } are symmetrically distributed j^. 

• By virtue of Theorems \l}and\2\ all methods in Examples U}(^nd\2\ are sym- 
metric, perfectly A-stable, and of order 2s. In particular such HBVM{k, s) 
are exact for polynomial Hamiltonian functions of degree v, provided that 

k>Y- (1-29) 
'That is, its region of Absolute stability precisely coincides with the left-half complex plane. 
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• For all k sufficiently large so that di.PD holds, HBVM{k, s) based on the 
k Gauss-Legendre abscissae in [0, 1] are equivalent to HBVM{k, s) based 
on k + 1 Lobatto abscissae in [0, 1] (see l\7]), since both methods define the 
same polynomial a of degree s ( i.e., they satisfy the same Master Functional 
Equation diJ5l) - 47?7?l) ). 



CHAPTER 1. BASIC FACTS ABOUT HBVMS 



Chapter 2 
Numerical Tests 



We here collect a few numerical tests, in order to put into evidence the potentiali- 
ties ofHBVMs flUillll. 

Test problem 1 

Let us consider the problem characterized by the polynomial Hamiltonian (4. 1) in 



having degree v = 6, starting at the initial point yo = (g(0),p(0))^ = (0, 1)^, so 
that Hiyo) = 0. For such a problem, in COl it has been experienced a numerical 
drift in the discrete Hamiltonian, when using the fourth-order Lobatto IIIA method 
with stepsize h = 0.16, as confirmed by the plot in Figure fTJl When using 
the fourth-order Gauss-Legendre method the drift disappears, even though the 
Hamiltonian is not exactly preserved along the discrete solution, as is confirmed 
by the plot in Figure [Z2l On the other hand, by using the fourth-order HBVM(6,2) 
with the same stepsize, the Hamiltonian turns out to be preserved up to machine 



precision, as shown in Figure [23l since such method exactly preserves polynomial 
Hamiltonians of degree up to 6. In such a case, according to the last item in 
Remark |71 the numerical solutions obtained by using the Lobatto nodes {cq = 
0, Ci, . . . , Ce = 1} or the Gauss-Legendre nodes {ci, . . . , Cg} are the same. The 
fourth-order convergence of the method is numerically verified by the results listed 



H{p, q) 




(2.1) 



in Table [2JJ 
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Figure 2.1: Fourth-order Lobatto IIIA method, h = 0.16, problem (|2.1|) : drift in 
the Hamiltonian. 



.X 10 




1500 



Figure 2.2: Fourth-order Gauss-Legendre method, h = 0.16, problem (12.11) : H 
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Figure 2.3: Fourth-order HBVM(6,2) method, h = 0.16, problem (lO) : H ^ 
10-16. 



17 



Test problem 2 

The second test problem, having a highly oscillating solution, is the Fermi-Pasta- 
Ulam problem (see [21, Section 1.5.1]), modelling a chain of 2m mass points con- 
nected with alternating soft nonlinear and stiff linear springs, and fixed at the end 
points. The variables gi, q2m stand for the displacements of the mass points, 
and Pi = Qi for their velocities. The corresponding Hamiltonian, representing the 
total energy, is 

H{p, 9) = 2 5Z (^2i-i + Pli)+^ '^^^i - <i2t~i f+Y^ (g2i+i - q2i)^ , (2.2) 

i=l i=l i=0 

with go = Q2m+i =0. In our simulation we have used the following values: 
m = 3, a; = 50, and starting vector 

Pi = 0, qi = {i-l)/lO, i = l,...,6. 

In such a case, the Hamiltonian function is a polynomial of degree 4, so that 
the fourth-order HBVM(4,2) method, either when using the Lobatto nodes or the 
Gauss-Legendre nodes, is able to exactly preserve the Hamiltonian, as confirmed 
by the plot in Figure IZ6l obtained with stepsize h = 0.05. Conversely, by using 
the same stepsize, both the fourth-order Lobatto IIIA and Gauss-Legendre meth- 
ods provide only an approximate conservation of the Hamiltonian, as shown in 
the plots in Figures and [231 respectively. The fourth-order convergence of the 
HBVM(4,2) method is numerically verified by the results listed in Table 12. 21 
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Figure 2.6: Fourth-order HBVM(4,2) method, h = 0.05, problem dO): \H - 
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Test problem 3 (non-polynomial Hamiltonian) 

In the previous examples, the Hamiltonian function was a polynomial. Neverthe- 
less, as observed above, also in this case HBVM(A;,s) are expected to produce a 
practical conservation of the energy when applied to systems defined by a non- 
polynomial Hamiltonian function that can be locally well approximated by a poly- 
nomial. As an example, we consider the motion of a charged particle in a magnetic 
field with Biot-Savart potential^ It is defined by the Hamiltonian [[61 

(2.3) 



with g = \/ -\- xp', a = e Bq, m is the particle mass, e is its charge, and Bq is 
the magnetic field intensity. We have used the values 

m = 1, e = —1, Bq = 1, 

with starting point 

a; = 0.5, y = 10, z = 0, x = -0.1, y = -0.3, i = 0. 

By using the fourth-order Lobatto IIIA method, with stepsize h = 0.1, a drift 
is again experienced in the numerical solution, as is shown in Figure 12.71 By 
using the fourth-order Gauss-Legendre method with the same stepsize, the drift 
disappears even though, as shown in Figure 12. 8[ the value of the Hamiltonian is 
preserved within an error of the order of 10^'^. On the other hand, when using 
the HBVM(6,2) method with the same stepsize, the error in the Hamiltonian de- 
creases to an order of 10^^^ (see Figure [Z91 ). thus giving a practical conservation. 
Finally, in Table 12.41 we list the maximum absolute difference between the numer- 
ical solutions over 10^ integration steps, computed by the HBVM(A;, 2) methods 
based on Lobatto abscissae and on Gauss-Legendre abscissae, as k grows, with 
stepsize h = 0.1. We observe that the difference tends to 0, as k increases. Fi- 
nally, also in this case, one verifies a fourth-order convergence, as the results listed 
in Table 12. 31 show. 



H{x,y,z,x,y,z) 



1 
2m 



X 



X — a- 



y~^~2) + + "log(^'))^ 



This kind of motion causes the well known phenomenon of aurora borealis. 
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Figure 2.7: Fourth-order Lobatto IIIA method, h = 0.1, problem (|2.3I) : drift in 
the Hamiltonian. 



5 



































ill 
















1 









500 1000 1500 2000 

t 

Figure 2.8: Fourth-order Gauss-Legendre method, h = 0.1, problem (12.31) : \H — 

^ xlQ-^^ 
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Figure 2.9: Fourth-order HBVM(6,2) method, h = 0.1, problem (1231) : \H - 
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Table 2.1: Numerical order of convergence for the HBVM(6,2) method, problem 



(1211). 



h 


0.32 0.16 0.08 0.04 0.02 


error 


2.288-10-2 1.487-10-3 9.398 - 10-^ 5.890 - 10-^ 3.684-10- 


7 


order 


3.94 3.98 4.00 4.00 


Table 2.2: Numerical order of convergence for the HBVM(4,2) method, problem 


(12.21). 






h 


1.6-10-2 8-10-3 4-10-3 2-10-3 10-3 




error 


3.030 1.967-10-1 1.240-10-2 7.761-10"^ 4.853 - 10"^ 




order 


3.97 3.99 4.00 4.00 




Table 2.3: Numerical order of convergence for the HBVM(6,2) method, problem 


(12.31). 






h 


3.2-10-2 1.6-10-2 8-10-3 4-10-3 2-10-3 


error 


3.944-10-6 2.635-10-^ 1.729-10"^ 1.094 - lO'^ 6.838-10" 


11 


order 


3.90 3.93 3.98 4.00 



Table 2.4: Maximum difference between the numerical solutions obtained through 
the fourth-order HBVM(A;, 2) methods based on Lobatto abscissae and Gauss- 
Legendre abscissae for increasing values of k, problem (12.31) . 103 steps with step- 
size /i = 0.1. 



k 


h = 0.1 


2 


3.97-10-1 


4 


2.29 - 10-3 


6 


2.01 - 10-s 


8 


1.37- 10-11 


10 


5.88 - 10-13 
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Test problem 4 (Sitnikov problem) 

The main problem in Celestial Mechanics is the so called A^-body problem, i.e. to 
describe the motion of N point particles of positive mass moving under Newton's 
law of gravitation when we know their positions and velocities at a given time. 
This problem is described by the Hamiltonian function: 

1 ^ ll« l|2 ^ ''^ m 

H{q, P) = ^E^-^E-^E (2.4) 

where is the position of the ith particle, with mass mj, and pi is its momentum. 

The Sitnikov problem is a particular configuration of the 3-body dynamics 
(see, e.g., OTTl ). In this problem two bodies of equal mass (primaries) revolve 
about their center of mass, here assumed at the origin, in elliptic orbits in the xy- 
plane. A third, and much smaller body (planetoid), is placed on the 2;-axis with 
initial velocity parallel to this axis as well. 

The third body is small enough that the two body dynamics of the primaries is 
not destroyed. Then, the motion of the third body will be restricted to the 2;-axis 
and oscillating around the origin but not necessarily periodic. In fact this problem 
has been shown to exhibit a chaotic behavior when the eccentricity of the orbits 
of the primaries exceeds a critical value that, for the data set we have used, is 
e ~ 0.725 (see Figure [2J0l) . 

We have solved the problem defined by the Hamiltonian function (12.41) by 
the Gauss method of order 4 (i.e., HBVM(2,2) at 2 Gaussian nodes) and by 
HBVM(18,2) at 18 Gaussian nodes (order 4, 2 fundamental and 16 silent stages), 
with the following set of parameters in (12. 4|) : 

G mi m2 e d h tmax 

3 111 10-5 0.75 5 0.5 1500 

where e is the eccentricity, d is the distance of the apocentres of the primaries 
(points at which the two bodies are the furthest), h is the used time-step, and 
[0, tmax] is the time integration interval. The eccentricity e and the distance d 
may be used to define the initial condition [qq, Pq] (see [|3T1l for the details): 

Qo = [-§, 0, 0, 1, 0, 0, 0, 0, 10-y, 

Po = [0, - j^VTo, 0, 0, ^VTo, 0, 0, 0, If. 

First of all, we consider the two pictures in Figure 12.1 II reporting the rela- 
tive errors in the Hamiltonian function and in the angular momentum evaluated 
along the numerical solutions computed by the two methods. We know that the 
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HBVM(18,2) precisely conserves Hamiltonian polynomial functions of degree at 
most 18. This accuracy is high enough to guarantee that the nonlinear Hamilto- 
nian function (12.41) is as well conserved up to the machine precision (see the upper 
picture): from a geometrical point of view this means that a local approximation 
of the level curves of (12.41) by a polynomial of degree 18 leads to a negligible er- 
ror. The Gauss method exhibits a certain error in the Hamiltonian function while, 
being this formula symplectic, it precisely conserves the angular momentum, as 
is confirmed by looking at the down picture of Figure 12.1 1[ The error in the 
numerical angular momentum associated with the HBVM(18,2) undergoes some 
bounded periodic-like oscillations. 

Figures 12.121 and 12.131 show the numerical solution computed by the Gauss 
method and HBVM(18,2), respectively. Since the methods leave the xy-plane 
invariant for the motion of the primaries and the z-axis invariant for the motion of 
the planetoid, we have just reported the motion of the primaries in the xy-phase 
plane (upper pictures) and the space-time diagram of the planetoid (down picture). 

We observe that, for the Gauss method, the orbits of the primaries are irregular 
in character so that the third body, after performing some oscillations around the 
origin, will eventually escape the system (see the down picture of Figure [2?T2l) . 
On the contrary (see the upper picture of Figure [2T3l) . the HBVM(18,2) method 
generates a quite regular phase portrait. Due to the large stepsize h used, a sham 
rotation of the xy-plane appears which, however, does not destroy the global sym- 
metry of the dynamics, as testified by the bounded oscillations of the planetoid 
(down picture of Figure 12.131) which look very similar to the reference ones in 
Figure l2?T0l This aspect is also confirmed by the pictures in Figure [OH display- 
ing the distance of the primaries as a function of the time. We see that the distance 
of the apocentres (corresponding to the maxima in the plots), as the two bodies 
wheel around the origin, are preserved by the HBVM(18,2) (down picture) while 
the same is not true for the Gauss method (upper picture). 
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Figure 2.10: The upper picture displays the configuration of 3-bodies in the Sit- 
nikov problem. To an eccentricity of the orbits of the primaries e = 0.75, there 
correspond bounded chaotic oscillations of the planetoid as is argued by looking 
at the space-time diagram in the down picture. 
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Figure 2.11: Upper picture: relative error \H(yn) — H(yQ)\/\H(yQ)\ oftheHamil- 
tonian function evaluated along the numerical solution of the HBVM(18,2) and 
the Gauss method. Down picture: relative error |M(y„) — M(yQ)\/\M{yo) \ of the 
angular momentum evaluated along the numerical solution of the HBVM(18,2) 
and the Gauss method. 
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Trajectory of 
Trajectory of m 





Figure 2.12: The Sitnikov problem solved by the Gauss method of order 4, with 
stepsize h = 0.5, in the time interval [0, 1500]. The trajectories of the primaries 
in the xy-plane (upper picture) exhibit a very irregular behavior which causes the 
planetoid to eventually escape the system, as illustrated by the space-time diagram 
in the down picture. 
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Figure 2.13: The Sitnikov problem solved by the HBVM(18,2) method (order 4), 
with stepsize h = 0.5, in the time interval [0, 1500]. Upper picture: the trajecto- 
ries of the primaries are ellipse shape. The discretization introduces a fictitious 
uniform rotation of the xy-plane which however does not alter the global symme- 
try of the system. Down picture: the space-time diagram of the planetoid on the 
z-axis displayed (for clearness) on the time interval [0, 350] shows that, although 
a large value of the stepsize h has been used, the overall behavior of the dynamics 
is well reproduced (compare with the down picture in Figure [2?T0l) . 
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Figure 2.14: Distance between the two primaries as a function of the time, related 
to the numerical solutions generated by the Gauss method (upper picture) and 
HBVM(18,2) (down picture). The maxima correspond to the distance of apoc- 
entres. These are conserved by HBVM(18,2) while the Gauss method introduces 
patchy oscillations that destroy the overall symmetry of the system. 



Chapter 3 
Infinity HBVMs 



From the previous arguments, it is clear that the orthogonality conditions (11.131) . 
i.e., the fulfillment of the Master Functional Equation (11.151) . is in principle only 
a sufficient condition for the conservation property (11.121) to hold, when a generic 
polynomial basis {Pj} is considered. Such a condition becomes also necessary, 
when such basis is orthonormal. 

Theorem 3. Let {Pj} be an orthonormal basis on the interval [0, 1]. Then, as- 
suming H{y) to be analytical, A1.12\) implies that each term in the sum has to 
vanish. 

Proof Let us consider the expansion 

g{T) = VH{a{to + rh)) = P^Pii^), Pi = (Pi, 9), ^ > 1, 

where, in general, 

if, 9)= [' f{r)g{r)dT. 
Jo 

Substituting into (11.121) . yields 

j=l j=l \ £>1 / j=l 

Since this has to hold whatever the choice of the function H{y), one concludes 
that 

7jPi = 0, j = l,...,s. □ (3.1) 
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Remark 8. In the case where {Pj} is an orthonormal basis, from d3.il) one then 
derives that 

-fj = Spj, i = l,...,s, 

where S is any nonsingular skew -symmetric matrix. The natural choice S = J 
then leads to A1.13\l . 

Moreover, we observe that, if the Hamikonian H{y) is a polynomial, the inte- 
gral appearing at the right-hand side in (11.181) is exactly computed by a quadrature 
formula, thus resulting into a HBVM(A;,s) method with a sufficient number of 
silent stages. As already stressed in the Chapter [H in the non-polynomial case 
such formulae represent the limit of the sequence HBVM(A;,s), as A; — oo. 

Definition 3. For general Hamiltonians, we call the limit formula A1.18\l Infinity 
Hamiltonian Boundary Value Method of degree s ( in short, oo-HBVM of degree 
s or HBVM(oo,s)j 

More precisely, due to the choice of the orthonormal basis (11.81) . 

HBVM(cx),s) = lim HBVM(A;,s), 

k—^oo 

whatever is the choice of the fundamental abscissae {q}. 

A worthwhile consequence of Theorems [T] and [2] is that one can transfer to 
HBVM(oo, s) all those properties of HBVM(A;,s) which are satisfied starting from 
a given k > ko on: for example, the order and stability properties. 

Corollary 1. Whatever the choice of the abscissae Ci, . . . , c^, HBVM{oo, s) (11.181) 
has order 2s and is perfectly A-stable. 



Chapter 4 



Isospectral Property of HB VMs and 
their connections with Runge-Kutta 
collocation methods 

When applied to initial value problems, HBVMs may be viewed as a special sub- 
class of Runge-Kutta (RK) methods of collocation type. In Chapter \T\ (see also 
(SI |7l) the RK formulation turned out useful in stating results pertaining to the 
order of the new formulae. Here, the RK notation will be exploited to derive the 
isospectral property of HBVMs and elucidate the existing connections between 
HBVMs and RK collocation methods [|9J. In doing this, our aim is twofold: 

1 . to better elucidate the close link between the new formulae and the classical 
collocation Runge-Kutta methods; 

2. to make the handling of the new formulae more comfortable to the scientific 
community working in the context of RK methods. 

In fact, we think that HBVMs (and consequently their RK formulation) may 
be of interest beyond their application to Hamiltonian systems. Each HBVM(A;,s) 
becomes a classical collocation method when k = s, while, for A; > s, it conserves 
all the features of the generating collocation formula, including the order (which 
may be even improved, reaching eventually order p = 2s) and the dimension of 
the associated nonlinear system. 

Let us then consider the matrix appearing in the Butcher tableau (11.271) . corre- 
sponding to HBVM(/c, s), i.e., the matrix 

A = IsVjn e M^■''^ k>s, (4.1) 

whose rank is s (see (|1.24l) - (l 1.251) ). Consequently it has a{k — s)-fold zero eigen- 
value. To begin with, we are going to discuss the location of the remaining s 
eigenvalues of that matrix. 
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Before that, we state the following preliminary result, whose proof can be 
found in ll23l Theorem 5.6 on page 83]. 



Lemma 1. The eigenvalues of the matrix 



6 







/ 



(4.2) 



with 







,s-l, 



(4.3) 



2v/(2j + l)(2j-l)' 

coincide with those of the matrix in the Butcher tableau of the Gauss-Legendre 
method of order 2s. 

We also need the following preliminary result, whose proof derives from the 
properties of shifted-Legendre polynomials (see, e.g., [1] or the Appendix in H). 



Lemma 2. With reference to the matrices in rti.24l) - di.25D . one has 

where 



(4.4) 



/ I -6 

6 



is-1 



(4.5) 



is J 



with the defined by t\4.3\) . 

The following result then holds true 

Theorem 4 (Isospectral Property of HBVMs). For allk > s and for any choice 
of the abscissae {rj} such that B(2s) holds true, the nonzero eigenvalues of the 
matrix A in M-.l^ coincide with those of the matrix of the Gauss-Legendre method 
of order 2s. 

Proof For k = s, the abscissae {rj} have to be the s Gauss-Legendre nodes 
on [0, 1], so that HBVM(s, s) reduces to the Gauss Legendre method of order 2s, 
as already observed in Example |2l 
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When k > s, from the orthonormality of the basis, see (11.71) . and considering 
that the quadrature with weights {cui} is exact for polynomials of degree (at least) 
2s — 1, one easily obtains that 



since, for all i 



s + 1: 



s, and 3 = 1,. 

By taking into account the result of Lemma [21 one then obtains: 

X, (/, 0) = n+iX, (/, 



AV 



s+1 



I 



\ -^1 












is-l 





is 


0/ 



v.. 



'Ps+lXs, 



(4.6) 



with the {^j} defined according to (14.31) . Consequently, one obtains that the 
columns of Vg+i constitute a basis of an invariant (right) subspace of matrix A, so 
that the eigenvalues of Xg are eigenvalues of A. In more detail, the eigenvalues of 
Xs are those of Xg (see (14.21) ) and the zero eigenvalue. Then, also in this case, the 
nonzero eigenvalues of A coincide with those of Xg, i.e., with the eigenvalues of 
the matrix defining the Gauss-Legendre method of order 2s. □ 



4.1 HBVMs and collocation methods 

By using the previous result and notations, now we go to elucidate the existing 
connections between HBVMs and RK collocation methods. We shall continue to 
use an orthonormal basis {-P,}, along which the underlying extended collocation 
polynomial a{t) is expanded, even though the arguments could be generalized to 
more general bases, as sketched below. On the other hand, the distribution of the 
internal abscissae can be arbitrary. 

Our starting point is a generic collocation method with k stages, defined by 
the tableau 



(4.7) 









A 


Tk 






UJl ... Uk 
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where, for i,j = 1, . . . ,k, A = (aij) = (/J"' £j(r)dr) andtuj = ij{T)dT, ij{t) 
being the jth Lagrange polynomial of degree k — 1 defined on the set of abscissae 

Given a positive integer s < A;, we can consider a basis {pi(r), . . . ,ps(r)} of 
the vector space of polynomials of degree at most s — 1, and we set 

/ Piin) P2{ri) ■■■ Psin) \ 



v.. 



(4.8) 



\ PliTk) P2iTk) ■■■ PsiTk) J ^ 



(note that Vs is full rank since the nodes are distinct). The class of RK methods 
we are interested in is defined by the tableau 



(4.9) 



Wl OJk 

where = diag(c<;i, . . . ,Uk) and = diag(?7i, . . . , 77^); the coefficients rjj, j = 
1, . . . ,s, have to be selected by imposing suitable consistency conditions on the 
stages {Yi} [TJ. In particular, when the basis is orthonormal, as we shall assume 
hereafter, then matrix Vs reduces to matrix Vs in (ll.24l) - (ll.25l) . = Ig, and 
consequently (14.91 ) becomes 



(4.10) 



n 






A = AVsVjn 


Tk 









We note that the Butcher array A has rank which cannot exceed s, because it 
is defined hy filtering A by the rank s matrix VsV^Vt. 

The following result then holds true, which clarifies the existing connections 
between classical RK collocation methods and HBVMs. 

Theorem 5. Provided that the quadrature formula defined by the weights {ui} is 
exact for polynomials at least 2s — 1 (i.e., the RK method defined by the tableau 
( \4.10\i satisfies the usual simplifying assumption B{2s)), then the tableau ( \4.10\) 
defines a HBVM{k, s) method based at the abscissae {rj}. 

Proof Let us expand the basis {Pi(r), . . . , Ps(r)} along the Lagrange basis 
{£j(r)}, j = 1, . . . , /c, defined over the nodes Ti, i = 1, . . . , k: 

k 

P,(r) = 5^P,(r,)4(r), J = l,...,s. 

r=l 
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It follows that, for i = 1, 



, k and j 




k 



k 



Pj{x)dx 



r=l 



r=l 



that is (see (fL24l) - (fL25l) and KT\ ). 



AVs. 



(4.11) 



By substituting (14.1 II) into (14.101) . one retrieves that tableau (|1.27l) . which defines 
the method HBVM(A;, s). This completes the proof. □ 

The resulting Runge-Kutta method (14.101) is then energy conserving if applied 
to polynomial Hamiltonian systems (|l.ll) when the degree of H(y), is lower than 
or equal to a quantity, say depending on k and s. As an example, when a 
Gaussian distribution of the nodes {rj} is considered, one obtains (11.291) . 

Remark 9 (About Simplecticity). The choice of the abscissae {ti, . . . ,Tk} at the 
Gaussian points in [0, 1] has also another important consequence, since, in such 
a case, the collocation method ( 14.71) is the Gauss method of order 2k which, as is 
well known, is a symplectic method. The result of Theorem\5\then states that, for 
any s <k, the HBVM{k, s) method is related to the Gauss method of order 2k by 
the relation: 



where the filtering matrix (VsVjQ) essentially makes the Gauss method of order 
2k "work" in a suitable subspace. 

It seems like the price paid to achieve such conservation properties consists in 
the lowering of the order of the new method with respect to the original one (14.71) . 
Actually this is not true, because a fair comparison would be to relate method 
(I1.27I) - (I4.10I) to a collocation method constructed on s rather than on k stages. 
This fact will be fully elucidated in Chapter [51 

4.1.1 An alternative proof for the order of HBVMs 

We conclude this chapter by observing that the order 2s of an HBVM(fc, s) method, 
under the hypothesis that (14.71) satisfies the usual simplifying assumption B(2s), 
i.e., the quadrature defined by the weights {coi} is exact for polynomials of degree 
at least 2s — 1, may be stated by using an alternative, though equivalent, procedure 
to that used in the proof of Theorem [T] 
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Let us then define the k x k matrix V = Vk (see (ll.24l) - (ll.25l) ) obtained by 
"enlarging" the matrix Vs with k — s columns defined by the normalized shifted 
Legendre polynomials Pjir), j = s + 1, . . . , k, evaluated at {tj}, i.e., 



V 



\ Pi{Tk) ... Pk{Tk) J 



By virtue of property B{2s) for the quadrature formula defined by the weights 
{coi}, it satisfies 



v^nv 



Is o 
O R 



Re 



i)k—sxk- 



This implies that V satisfies the property T(s, s) in [|23l Definition 5. 10 on page 
86], for the quadrature formula (ui, Tj)^^]^. Therefore, for the matrix A appearing 
in (14.101) (i.e., (11.271) . by virtue of Theorem [5]), one obtains: 



V-^AV = V-^AV 



O 



o 



(4.12) 



where Xg is the matrix defined in (14.61) . Relation (14.121) and J23l Theorem 5.1 1 on 
page 86] prove that method (14.101) (i.e., HBVM(A;, s)) satisfies C{s) and D{s - 1) 
and, hence, its order is 2s. 

Remark 10 (Invariance of the order). From the previous result we deduce the 
invariance of the superconvergence property of HBVM(k,s) with respect to the 
distribution of the abscissae Ti, i = 1, . . . , k, the only assumption to get the order 
2s being that the underlying quadrature formula has degree of precision 2s — 1. 
Such exceptional circumstance is likely to have interesting applications beyond 
the purposes here presented. 



Chapter 5 
Blended HBVMs 



We shall now consider some computational aspects concerning HBVM(/c, s). In 
more details, we now show how its cost depends essentially on s, rather than on 
k, in the sense that the nonlinear system to be solved, for obtaining the discrete 
solution, has (block) dimension s ["SlISllHl- 

This could be inferred from the fact that the silent stages (|1.1 II) depend on the 
fundamental stages: let us see the details. In order to simplify the notation, we 
shall fix the fundamental stages at ri, . . . , r^, since we have already seen that, due 
to the use of an orthonormal basis, they could be in principle chosen arbitrarily, 
among the abscissae {rj}. With this premise, we have, from (11.91) . (I1.17I) - (I1.18I) . 
and by using the notation (11.211) . 

s k 

Vi = yo + h^a.i,j^uJiPj{Te)fi, i = l,...,s. (5.1) 

j=i e=i 

This equation is now coupled with that defining the silent stages, i.e., from 
(fL6l) and (fLTTT) . 

yi = yo + hy^'yj / Wd^) i = s + l,...,k. (5.2) 
Let us now partition the matrices Is,Vs E R''^'' in (ll.24l) - (ll.25l) into 

containing the entries defined by the fundamental abscissae and the silent ab- 
scissae, respectively. Similarly, we partition the vector y into yi, containing the 
fundamental stages, and ^2 containing the silent stages and, accordingly, let 

Qi E M'^^ ^2 E Ri'-'x^-^ 
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be the diagonal matrices containing the corresponding entries in matrix fi. Finally, 
let us define the vectors 

7 = (71, . . . , 7,)^, e = (1, . . . , 1)^ G M^ u={l,..., if e R'-'. 

Consequently, we can rewrite (15.11) and (15.21) . as 

= e®yo + hXs, [Vl Vj,) (^^' ^ J ® hm ( If^^l ) , (5.3) 

1/2 = u®yQ + hXs2® hml, (5.4) 
respectively. The vector 7 can be obtained by the identity (see (|1.16|) ) 

t/i = e (g) yo + hXsi ® hml, 

thus giving 

= (m - X,2X7/e) ® 1/0 + ^,2X7/ (g) hmVi 

= u^yo + Aii^ hmVi, (5.5) 

in place of (15.41) . where, evidently, 

u={u- Xs2X:^e) e M'=-^ = X,2X7/ G R'^-^^^ (5.6) 

By setting 

B, = XsiVlQi e R^^^ B2 = XsiVj^Q2 e M^^'=-^ (5.7) 



substitution of (15.51) into (15.31) then provides, at last, the system of block size s to 
be actually solved: 

F{yi) = yi-e^yo-h[Bi^l2mf{yi)+ (5.8) 

B2 ® hmf (m ® 2/0 + ^1 ® hmVl)] = 0. 

By using the simplified Newton method for solving (15.81) . and setting 

C = Bi + B2A1 G M'^^ (5.9) 

one obtains the iteration: 

{h®l2m-hC®J,)5^^^ = -F{y^r^) = ^p^r\ (5.10) 



(n+l) 

Vi = y 



S") + 5("), n = 0,l,... 



where Jo is the Jacobian of f{y) evaluated at yo. Because of the result of Theo- 
rem IH the following property of matrix C holds true |[8l. 
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Theorem 6. The eigenvalues of matrix C in f l5.9l) coincide with those of matrix 
( \4.2\i . i.e., with the eigenvalues of the matrix of the Butcher array of the Gauss- 
Legendre method of order 2s. 

Proof Assuming, as usual for simplicity, that the fundamental stages are the 
first s ones, one has that the discrete problem 

^ " ( M ) ^Vo + hA^h^fiy), 

which defines the Runge-Kutta formulation of the method, is equivalent, by virtue 
of (O, dSJ), dUl), dH]), to 

Is Ogxr \ ^ T I Vl 



where, as usual, r = k — s. Consequently, the eigenvalues of the matrix A defined 
in (|4.1I) coincides with those of the pencil 



That is. 



~Ai If J \ Orxs Oj-xi 



-Ai Ir J \ V J I Orxs Orxr ) \ V 



(5.11) 



for some nonzero vector (w^, -y^)^. By setting u = 0, one obtains the r zero 
eigenvalues of the pencil. For the remaining s (nonzero) ones, it must hev = Aiu, 
so that: 

Hu= {Biu + B2V) = {Biu + B2A1U) = Cu ^ G o-(C). □ 



Remark 11. From the result of Theorem^ it follows that the spectrum of C 
doesn 't depend on the choice of the s fundamental abscissae, within the nodes 
{ti}. On the contrary, its condition number does: the latter appears to be mini- 
mized when the fundamental abscissae are symmetrically distributed and approx- 
imately evenly spaced in the interval [0, 1]. As a practical "rule of thumb", the 
following algorithm appears to be almost optimal: 

1. let the k abscissae {rj} be chosen according to a Gauss-Legendre distribu- 
tion ofk nodes; 
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2. then, let us consider s equidistributed nodes in (0, 1), say {fi, . . . , f^}; 

3. select, as the fundamental abscissae, those nodes among the {tj} which are 
the closest ones to the {r-, }; 

4. define matrix C in d5. 91) accordingly. 

Clearly, for the above algorithm to provide a unique solution ( resulting in a sym- 
metric choice of the fundamental abscissae), the difference k — s has to be even 
which, however, can be easily accomplished. 

In order to give evidence of the effectiveness of the above algorithm, in Fig- 
ure [5TT] we plot the condition number of matrix C = C{k,s), for s = 2, ... ,5, 
and A; > s. As one can see, the condition number of C (k, s) turns out to be nicely 
bounded, for increasing values of k, which makes the implementation (that we are 
going to analyze in the next section) effective also when finite precision arithmetic 
is used. For comparison, in Figure 15.21 there is the same plot, obtained by fixing 
the fundamental abscissae as the first s ones. In such a case, the condition number 
of C{k, s) grows very fast, as k is increased. 

5.1 Blended implementation 

We observe that, since C is nonsingular, we can recast problem (15.101) in the equiv- 
alent form 

7 {C-' ® hn. - hi, ® Jo) S^^'^ = ® hrnFiy^^) = (5.12) 

where 7 > is a free parameter to be chosen later. Let us now introduce the 
weight (matrix) function 

d = h®^-\ $ = /2„ - /i7Jo G M""'''"', (5.13) 
and the blended formulation of the system to be solved, 

M5(") = [d{Is®hm-hC ® Jq) + 

(/ - ^)7 (C-l ® hm - his ® Jo)] S^''^ 

= ^vJ"^ + {1- = V'^"^- (5.14) 

The latter system has again the same solution as the previous ones, since it is 
obtained as the blending, with weights and (1 — 9), of the two equivalent forms 
(15.101) and (15.121) . For iteratively solving (15.141) . we use the corresponding blended 
iteration, formally given by|Il[ia[IIl[I3[Il[l4l[l5l[l6l[l8l[30l[3l: 

^(M+i) = ^(M) _ ^ (^jVf^M _ j ^ £ = 0,1,.... (5.15) 



5. 1 . BLENDED IMPLEMENTATION 
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Figure 5.2: Condition number of the matrix C = C{k, s), for s = 2, 3, 4, 5 and 
k = s, s + 1, . . . , 100, with the fundamental abscissae chosen as the first s ones. 
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Remark 12. A nonlinear variant of the iteration d5.i5D can be obtained, by start- 
ing at S^'"''^'' = and updating il^^"^^ as soon as a new approximation is available. 
This results in the following iteration: 

^(n+l) _ ^(n) ^^^(n)^ ^ = 0,1,.... (5.16) 

Remark 13. We observe that, for actually performing the iteration ( \5.13\i - ^5.15\i . 
as well as 0.7(51) . one has to factor only the matrix $ in l\5.13\) . which has the same 
size as that of the continuous problem. 

We end this section by observing that the above iterations (15.151) and (15.161) 
depend on a free parameter 7. It will be chosen in order to optimize the conver- 
gence properties of the iteration, according to a linear analysis of convergence, 
which is sketched in the next section. 



5.2 Linear analysis of convergence 

The linear analysis of convergence for the iteration (15.151) is carried out by con- 
sidering the usual scalar test equation (see, e.g., ^141 and the references therein), 

y' = \y, 3?(A) < 0. 

By setting, as usual q = hX, the two equivalent formulations (15.101) and (15.121) 
become, respectively (omitting, for sake of brevity, the upper index n), 

{Is - qC)6 = 7(C-i - qls)d = li,^. 

Moreover, 

e = e{q) = {l-jqr'h, (5.17) 
and the blended iteration (15.151) becomes 

5(^+1) = (/^ _ 9{q)M{q))d^'^ + 9{q),/,{q), (5.18) 

with 

M(g) = 9iq) - qC) + {Is ~ 9iq))j - qls) , (5.19) 

V'(g) = ^(g)V'i + (/s - ^(g))V2- 

Consequently, the iteration will be convergent if and only if the spectral radius 
p{q) of the iteration matrix, 

Z{q)=h-e{q)M{q), (5.20) 

is less than 1. The set 

T = {qeC: p{q) < 1} 
is the region of convergence of the iteration. The iteration is said to be: 
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Table 5.1: optimal values (15.231) . and corresponding maximum amplification fac- 
tors (15.241) . for various values of s. 



s 


7 


P* 


2 


0.2887 


0.1340 


3 


0.1967 


0.2765 


4 


0.1475 


0.3793 


5 


0.1173 


0.4544 


6 


0.0971 


0.5114 


7 


0.0827 


0.5561 


8 


0.0718 


0.5921 


9 


0.0635 


0.6218 


10 


0.0568 


0.6467 



• ^-convergent, if C C F; 

• L-convergent, if it is A-convergent and, moreover, p(g) — )■ 0, as g — )• oo. 
For the iteration (15.181) one verifies that (see (15.171) . (15.191) . and (|5.20l) ) 

Z{q) = ^^^.C-' {C - 7/.)' , (5.21) 

which is the null matrix at g = and at oo. Consequently, the iteration will be A- 
convergent (and, therefore, L-convergent), provided that maximum amplification 
factor, 

p* = max p(q) < 1. (5.22) 

SR(g)=0 

From (|5.21l) one has that, by setting hereafter cr(C) the spectrum of matrix C, 
By taking into account that 

\q\ 1 

max 



SR(5)=o 1(1 -7g)2| 27' 
one then obtains that 



p = max 



\P-1\ 



For Gauss-Legendre methods (and, then, for any matrix C having the same spec- 
trum), it can be shown that (see [fT0l[T6ll ) the choice 

7 = \Pmm\ = min (5.23) 

tJ.&cr{C) 
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minimizes p*, which turns out to be given by 

p* = 1- COSV^min < 1, V^min = Arg(/imm)- (5.24) 

In Table I5.1[ we list the optimal value of the parameter 7, along with the 
corresponding maximum amplification factor p*, for various values of s, which 
confirm that the iteration (|5.18l) is L-convergent. 

Remark 14. We then conclude that the blended iteration 0.751) turns out to be 
L-convergent, for any HBVM{k, s) method, for all s > 1 and k > s. 

We end this chapter, by emphasizing that the property of L-convergence has 
proved to be computationally very effective, as testified by the successful imple- 
mentation of the codes BiM and BiMD [[30l l32l . We then expect good perfor- 
mances also for the blended implementation of HBVM(A;, s). 



Chapter 6 

Notes and References 



The approach of using discrete line integrals has been used, at first, by lavernaro 
and Trigiante, in connection with the study of the properties of the trapezoidal rule 

mmm. 

It has been then extended by lavernaro and Pace ll24l . thus providing the first 
example of conservative methods, basically an extension of the trapezoidal rule, 
named s-stage trapezoidal methods: this is a family of energy-preserving methods 
of order 2, able to preserve polynomial Hamiltonian functions of arbitrarily high 
degree. 

Later generalizations allowed lavernaro and Pace [|25l . and then lavernaro and 
Trigiante [29|, to derive energy preserving methods of higher order. 

The general approach, involving the shifted Legendre polynomial basis, which 
has allowed a full complete analysis of HBVMs, has been introduced in (611 (see 
also jSl) and, subsequently, developed in [|3. 

The Runge-Kutta formulation of HBVMs, along with their connections with 
collocation methods, has been studied in fOl. 

The isospectral property of HBVMs has been also studied in [8J, where the 
blended implementation of the methods has been also introduced. 

Computational aspects, concerning both the computational cost and the effi- 
cient numerical implementation of HBVMs, have been studied in |[3l and [8]. 

Relevant examples have been collected in [4J, where the potentialities of HB- 
VMs are clearly outlined, also demonstrating their effectiveness with respect to 
standard symmetric and symplectic methods. 

Blended implicit methods have been studied in a series of papers ^ [TOl dH 
[121 [131 [Ml [151 [161 [30l and have been implemented in the two computational codes 
BiM and BiMD [l32l . 
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